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Obtaining a rigorous and reliable method for linking computer simulations of polymer blends and 
composites at different length scales of interest is a highly desirable goal in soft matter physics. 
In this paper a multiscale modeling procedure is presented for the efficient calculation of the static 
structural properties of binary homopolymer blends. The procedure combines computer simulations 
of polymer chains on two different length scales, using a united atom representation for the finer 
structure and a highly coarse-grained approach on the meso-scale, where chains are represented 
as soft colloidal particles interacting through an effective potential. A method for combining the 
structural information by inverse mapping is discussed, allowing for the efficient calculation of 
partial correlation functions, which are compared with results from full united atom simulations. 
The structure of several polymer mixtures is obtained in an efficient manner for several mixtures 
in the homogeneous region of the phase diagram. The method is then extended to incorporate 
thermal fluctuations through an effective \ parameter. Since the approach is analytical, it is fully 
transferable to numerous systems. 



I. INTRODUCTION 

Mixtures of polymers with different compositions are 
of great interest for many technological and industrial ap- 
plications. For example, multicomponent polymer mix- 
tures facilitate the custom tailoring of desired physical 
properties, leading to the creation of new materials with 
enhanced performance. pQ However, since new materials 
are typically processed in their liquid phase, the ability 
to predict the mixing behavior of various polymer com- 
posites is a desirable computational task. Because of the 
underlying phase diagram, polymer mixtures can phase 
separate at specific thermodynamic conditions of temper- 
ature and polymer volume fraction. The structure and 
dynamics of each component strongly depends on the re- 
gion of the phase diagram that is explored, and different 
length scales characterize the thermodynamic and dy- 
namical behavior of the mixture, depending on how far 
the system is from its phase separation. 

While computer simulations are useful in elucidating 
the structural and dynamic complexity of these systems, 
they are limited in the extent of information that can 
be collected by the precision of the calculation, which 
degrades with the number of simulation steps, and by 
the power of the available computational machines. [2] Al- 
ready for polymer liquids it is not possible to collect all 
the relevant information from a single simulation when 
long polymers are involved, as information spreads over 
many order of magnitude in time and length, going from 
the monomer length scale, I, to the overall polymer di- 
mension of the radius-of-gyration, R g = y/Nl/6. The 
latter scales as the square root of the degree of poly- 
merization, N, which can be fairly large for polymers of 
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interest. [3] The added complexity of dealing with poly- 
mer mixtures arises from the presence of a new length 
scale of concentration fluctuations, which diverges 
as the system approaches its spinodal curve. The big 
challenge for such systems becomes simulating the mix- 
tures not only at different temperatures and composi- 
tions, but also with a box size that increases with the 
scale of concentration fluctuations. Furthermore, the 
equilibration step in molecular dynamics (MD) simula- 
tions becomes extremely slow, as the time scale in which 
polymer mixtures phase separate diverges as the temper- 
ature approaches the critical point. Finally the timescale 
of relaxation of the mixture depends on the proximity of 
the system to the glass transition of its components. In 
conclusion, classical MD simulations of polymer mixtures 
have been limited so far to high temperatures close to the 
athermal limit. In that region the thermodynamic be- 
havior is mainly guided by entropy and local packing[3], 
while enthalpic effects may be accounted for within a high 
temperature perturbative framework. [?J [5] 

In order to extend the ability of MD simulations to 
make quantitative predictions about polymer blends at 
temperatures approaching the spinodal, information is 
needed on multiple length-scales of interest. As the 
system approaches its phase transition and the largest 
length scale of interest diverges, one might be tempted 
to neglect the details of the local scale; however, such 
a procedure risks losing pertinent information, since the 
monomeric structure of a blend's component largely de- 
termines its glass transition temperature, as well as the 
shape of its phase diagram. Nonetheless, capturing the 
global structure is equally as important since an accurate 
determination of the thermodynamics of the system re- 
quires obtaining the total distribution of particles in the 
ensemble. For example, while system-specific local infor- 
mation is contained within the correlation hole region of 
the pair distribution function, g(r), the global structure 
is still required, as truncation of g(r) will prevent accu- 
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rate Fourier transform to the structure factor where the 
low k limit is important. (6] Thus, both local and global 
scale information need to be accounted for to achieve a 
complete thermodynamic and dynamic description of the 
mixture. 

One strategy to improve our ability to simulate mix- 
tures of complex fluids and collect information on an in- 
creasingly large range of length scales is to use a hier- 
archical approach, such as the multiscale modeling pro- 
cedure that we present in this paper. In such a multi- 
scaling scheme the system of interest is represented at 
various levels of coarse-graining and simulations are per- 
formed for each of these coarse-grained representations to 
enhance computational efficiency. In a second step, the 
information obtained from those simulations is collected 
and combined into a complete description, bridging in- 
formation at all the length scales of interest. [7] 

For any coarse-grained model to be capable of faith- 
fully reproducing the desired characteristics of the sys- 
tem, the coarse level of description needs to include the 
essential features of the atomistic structure while aver- 
aging out the remaining degrees of freedom. Several 
coarse-graining procedures have been developed in the 
literature at various levels of "coarseness" . At the lowest 
level, the united atom (UA) description, which represents 
groups of CH X units, with x = 1, 2, or 3 into an effective 
UA along the chain, is a particularly useful represen- 
tation, for which potentials have been optimized to re- 
produce the structure for most polyolefins.[8 Higher lev- 
els of structure-based approaches simplify interactions, 
while remaining close to the chemical structure, facilitat- 
ing scale hopping between the different levelsJS]; however, 
these methods remain computationally demanding and 
do not adopt an extreme enough level of coarse-graining 
to mitigate the limitations of fully atomistic simulations 
in reaching the large length scales of phase separation. [TU] 
In order to simulate the increasingly large scale phenom- 
ena of interest to material scientists, an even coarser level 
of description is needed, while remaining formally related 
to the finer structure. 

The highest level of coarse-graining in our procedure, 
represents each polymer chain as a soft colloidal particle, 
centered at the polymer center-of-mass. Soft colloids in- 
teract through a Gaussian-like repulsive potential of the 
order of the polymer radius-of-gyration, R g , which de- 
fines the characteristic length scale of this description. 
Large-scale behavior is easily generated in a "mesoscale" 
simulation of the coarse-grained, mixture and informa- 
tion is collected on a scale equal to or larger than the 
polymer radius-of-gyration. Recently, we have shown 
that our mesoscale simulations capture the global struc- 
ture of the polymer mixture providing information over 
much of the miscible region of the phase diagram. [llj 

Much attention has been given recently to relating 
mesoscopic [P2~j or field-theoretical model parameters [13] 
to their more detailed atomistic representation. For ex- 
ample, Groot and Warren were able to express the re- 
pulsive soft potential parameters used in dissipative par- 



ticle dynamics (DPD) in terms of the Flory-Huggins \ 
parameter, providing a bridge between DPD simulations 
and atomistic simulations, from which solubility parame- 
ters can be directly calculated. [2] This method has been 
successfully implemented in mesoscale simulations of bi- 
ological membranes once the model is parameterized for 
specific molecular structures. [15] In our approach, the ef- 
fective x parameter is included directly in the model as an 
input parameter which determines the length scale of the 
of concentration fluctuations. This \ parameter provides 
a convenient link between an atomistic or united atom 
level description and a mesoscale picture, similar to other 
mesoscale modeling approaches: [TBI ITU ITC] however, our 
approach has the advantage over other methods such as 
DPD in that the soft-core potential depends directly on 
the x parameter and not on any phenomenological ex- 
pression to bridge the two theories. In other words, once 
the x parameter is specified, our multiscale modeling pro- 
cedure will produce the expected phase behavior [TTJ, and 
the task in practice is to determine this parameter for a 
given model. 

In the current paper, we focus on bridging the length 
scales of binary homopolymer mixtures of different chem- 
ical architecture to provide a complete description of the 
static structure of the blend under a given set of ther- 
modynamic conditions. Blends of polyolefins represent a 
good test case for our procedure since united atom pa- 
rameters have been well optimized for these systems, and 
they provide a stringent test of the method, as subtle 
differences in chemical architecture, such as the extent of 
branching, greatly affects the global structure and mis- 
cibility. One advantage of our coarse-graining and mul- 
tiscale modeling approaches is that since they are based 
on the formalism of the Ornstein-Zernike equation, the 
rescaling of the structure to recover the monomer level 
of description is straightforward. Furthermore, we have 
shown that the correct dynamics can be obtained di- 
rectly from simulations of the coarse-grained structures 
by proper rescaling. |17j Thus, by performing simulations 
of mesoscopic particles and subsequently reinserting the 
relevant chemical details afterwards, we are able to ob- 
tain quantitative information about large-scale proper- 
ties. 

Large scale information, obtained form mesoscale MD 
simulations, is later combined with local-scale united- 
atom simulations of the same system, using a multiscale 
formalism obtained from the solution of the Ornstein- 
Zernike equation, presented in this paper. In a spirit 
similar to the united atom simulations to which we test 
our approach, we initially take the high temperature limit 
as a reference state, and set x = in order to calculate 
monomer level total correlation functions. The model is 
then extended to incorporate thermal effects on the con- 
centration fluctuation through the interaction parameter, 
X, which serves to bridge the gap between mesoscale and 
atomistic models. 

This multiscale modeling approach, is an extension 
to polymer mixtures of our coarse-graining and multi- 



3 



scale modeling procedure for polymer melts. [TS] By com- 
parison with united atom simulations of polymer mix- 
tures, our approach shows good agreement between the 
two descriptions. Furthermore, our procedure has the 
advantage of being analytically solved, so that thermo- 
dynamic parameters defining the free energy potential 
appear explicitly in the formalism, making the whole 
procedure general and transferible to other systems in 
different thermodynamic conditions. The key advan- 
tage of our method lies in the fact that it provides a 
quantitative picture of how the local details such as the 
specific chemical architecture and thermodynamic con- 
ditions are connected to the observed, global properties. 
While atomistic and united atom simulations are limited 
in the maximum length scale that can be simulated, our 
approach easily provides information on a wide range of 
length scales, once mesoscale simulations are combined 
with atomistic or monomer level simulations. The num- 
ber of test systems presented in the paper is determined 
by the available united-atom simulations against which 
our procedure is tested; however, our procedure is di- 
rectly applicable to other mixtures as well. 

While the focus here is on the equilibrium structural 
properties of simple polyolefins, it is worth noting that 
many interesting systems could be investigated through 
our method and important problems approached. For 
example, the proposed multiscale procedure could be ap- 
plied to study various aspects of binary blends, including 
critical behavior, mechanical properties and viscosity, the 
effects of non-random mixing, stiffness disparity, shape 
disparity, and different architectures, as well as the dy- 
namics of demixing at equilibrium and under shear, in- 
tcrdiffusion in equilibrium and during phase separation, 
and the different propensity for crystallization in the dy- 
namics of miscible blends approaching the phase transi- 
tion. As one of the polymers in the mixture crystalizes 
it is important to investigate, for example, the kinetics 
of interface formation and structure. Additionally the 
procedure could be extended to encompass mixtures in- 
cluding block and homo-copolymers, blends close to sur- 
faces and in thin films, and mixtures of colloids and poly- 
mers. Although the dynamics of coarse-grained systems 
is unrealistically accelerated, we have recently proposed 
a new first-principle method to rescale the dynamics of 
the coarse-grained system to its atomistic value, enabling 
dynamical properties of these systems to be investigated 
with quantitative precision. [17] This method can be use- 
ful in evaluating long-time dynamics directly from the 
mesoscale simulation. 

The paper begins with a brief overview of the the- 
ory for coarse-graining of polymer blends, from which 
we calculate the effective potential, v(r), input to the 
mesocale molecular dynamic simulations of the coarse- 
grained system. Mesoscale simulations of the coarse- 
grained polymer mixture, where chains are represented 
as soft colloidal particles, are presented in the following 
section. Those simulations provide the large-scale infor- 
mation, which is combined with local-scale united-atom 



simulations by means of our multiscale modeling proce- 
dure. Next, we present our results in the form of the 
total correlation function, h a< p(r) with a, @ £ A, B for 
each pair of components, which is predicted and com- 
pared with data of full united-atom simulations. First, 
we present correlation functions of the mixture in the 
limit where x = to demonstrate the tenor of our scale 
bridging approach. This is followed by an extension of 
the model to include finite temperature effects, where we 
present the change in the concentration structure factor 
upon inclusion of realistic values for the \ parameter. A 
brief discussion concludes the paper. 



II. METHODOLOGY 

A. Effective Pair Potentials for Polymer Mixtures 

An analytical coarse graining procedure for binary poly- 
mer blends was derived by Yatsenko, et aL[19], to for- 
mally map the blend onto a mixture of soft-colloidal 
particles. Here, we briefly review the results of that 
work. The polymer species (A and B) in the blend are 
characterized by a radius of gyration, R 9 a and R 9 b, 
a total density of monomers, p, and a segment length 
oa = \/Q/NARgA and <jb = y/^/NBRgB, where Na and 
Nb are the number of chain units. The volume fraction 
of the A component is given by <f> and the ratio of seg- 
ment lengths is given by 7 ~ ga/^b- The main physical 
quantity of interest in the study of complex liquids, is the 
pair distribution function, or the related total correlation 
function, as from this quantity any thermodynamic prop- 
erty of interest can be directly calculated. [5U] Here the 
total correlation function has contributions from three 
terms: hAA(k), /iss(fc), fiAB(k), representing the three 
types on interaction of a mixture of particles A and B. 
Analytical expressions for the total correlation functions 
at the coarse grained level are obtained by solving a gen- 
eralized Ornstcin-Zernike (OZ) matrix relation 



HO) = n(k)C(k) [fi(fc) + H(fc)] 



(1) 



where H(k) is the matrix of total intermolecular partial 
correlation functions (pefs), C(k) is the matrix of direct 
pefs, and represents the matrix of intramolecular 

pefs. Solving Equation [T] gives the relation for a generic 
pair[H] {a,p & (A, B)} 



h%{k) 



r-mm 
tl a(3 



(k), 



(2) 



where uj{k) for each type of chain is the intramolecular 
form factor, and the superscripts, cc, cm, or mm, de- 
notes center of mass-center of mass, monomer-center of 
mass, or monomer-mononomer interactions, respectively. 
The derivation of Equation [2] is an extension of the pro- 
cedure outlined by Krakoviack, et al.\2.2\ for homopoly- 
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mers and formally connects the center of mass (cm) dis- 
tribution functions to monomer-monomer (mm) distri- 
bution functions. Both the intramolecular form factors 
in Equation [2j namely the monomer- monomer uj mm {k) 
and the monomer-cm io cm (k), are expressed in analyt- 
ical forms. [23] [24] Yatsenko, et al. introduced a new 
analytical form for the monomer total correlation func- 
tions, h™™ by extending the PRISM-blend thread model 
of Tang and Schweizer|25) to include asymmetries in local 
chemical structure and flexibility. 19 

Following this procedure, analytical forms of h^ s , 
coarse-grained at the cm level, can be calculated from 
Equation [2] In real space they are given by: 

1 



(r) 



-lUr)+l 2 I P AA {r) 



h B B{r) 
h'AB (r) 



^M + ^W, (3) 
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parameter which determines the length scale of concen- 
tration fluctuations. The choice of a single adjustable 
parameter is appropriate in this case, since the aim is to 
obtain quantitative agreement with experiment, and the 
adoption of a single "effective" x maintains a straight- 
forward connection to SANS experiments. The notion of 
extracting a single effective x parameter was also invoked 
Dudowicz et al. in a general lattice cluster theory anal- 
ysis of compressible blends, [57] and by Schweizer in con- 
necting the PRISM blend approach with methods used 
in experimental SANS analysis. [35J While we initially set 
X = 0, the extension to more realistic models is achieved 
by using the experimental x parameter within this gen- 
eral framework, which has been shown to represent well 
both upper and lower critical phase diagrams. |21) Results 
presented in Section III B of this paper are examples of 
both kinds of systems. In the athermal regime, however, 
the theory is presently being implemented to reproduce 
quantitatively and in a self-consistent way, the experi- 
mental values of Xath [2S] ■ 

Finally, it should be noted that here we are only inter- 
ested in capturing the global structure and thus the use 
of the simple analytic PRISM thread result is adequate 
as it accurately captures the structure of the polymer at 
(4)distances greater the R g and correctly predicts the cor- 
relation hole. While the PRISM thread result does not 
capture local effects such as solvation shells, such detail is 
not relevant at the coarse-grained level since at this level, 
local interactions are averaged out. Local packing in the 
structural g(r) is introduced later using the multiscale 
modeling procedure discussed later, where local scale in- 
formation from UA MD simulations is paired with global 
information obtained from course-grained simulations. 

By adopting the hypernetted-chain (HNC) closure re- 
lation, the effective pair interaction potential, V^ g (r), is 
connected to the pefs by the relationship 



where £, ca p is the average correlation hole length scale, (^bT) 1 v a c g (r) — h L a c g (r) ln[l + h^ g (r)] c a c 8 (r) , (8) 

where c c a c g is the direct pcf for a mixture of soft colloidal 
particles, defined in reciprocal space as 



and Rga/3 = ^ {R 2 ga +R 2 g/3 )/2 = CcapV^ is an average 

radius-of-gyration. Here, £a € {£p,£$} identifies the 
length scale for the density and concentration fluctuation 
correlations. The concentration fluctuation of a polymer 
blend is defined as 

CAB 



(7) 



1 

Pea (Pc,a + Pc,/3)|S cc (fc)| 
(Pc,a + Pc,p)\Scc{fy\ 



(9) 



where we introduce the miscibility parameter x for which 
the quantity x/p is analogous to the Flory-Huggins in- 
teraction parameter [21] such that x The con- 
centration fluctuation, diverges at the spinodal tem- 
perature where x — Xs- The average segment length is 
a 2 AB = 4><r B + (l — <p)(y\, while the density fluctuation, £ p , 
is defined as ^ = npal B /3+^ and £' p = Z / (j: pa\ B ). 
Although for compressible blends, there is, in princi- 
ple, more than one x parameter resulting from the dis- 
tinct partial osmotic compressibilities of a two compo- 
nent blend [55]. in this work we adopt a single interaction 



The chain density, p c , a , of chain type, A, is given as 
Pc,A = 4>p/Na, and for chain type, B, as p c ^ B = (1 — 
4>)p/Nb- For a binary mixture the static structure fac- 
tors, S a c g , are given by 



Saa(Q = <t> + 4> 2 pch h A A(k) , 
S B C B(k) = l-^+(l-(b) 2 p ch h BB (k) 

SAB(k) = 0(l-$PcfcW*) . 



(10) 
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where the total chain density p c h — p/N, and |S cc (fc)| = 
S A c A (k)S% B {k) - [S c A c B (k)} 2 is the determinant of the 
mesoscopic static structure factor matrix. Whereas at 
the monomer level of description, molecular closures 
are required to ensure the correct scaling of the % 
parameter |30j. the use of the site level HNC is appropri- 
ate here because at this level of coarse-graining, the poly- 
mers are treated as a simple liquid of soft-colloids. 31 

By inserting Equations [3] and [9] into equation [8j the 
effective pair potentials are calculated, which are input 
to the mesocale simulations of the coarse grained binary 
mixture. The potential so obtained explicitly depends 
on the structural parameters of the polymer, such as p, 
N, cj>, and R g . Since these parameters enter into the 
UA description, they do not represent additional pa- 
rameters needed in moving to a higher level of coarse- 
graining. In other words, the potential is state depen- 
dent, being a free energy obtained from the monomer 
frame of reference; however, it is fully transferable to 
different systems since it is analytically derived in a well- 
defined manner. In addition to these structural param- 
eters, there is the additional parameter, x> that enters 
into the mesoscale description through Equation [7j and 
describes the monomer-specific interactions that drive 
phase separation. Once these parameters are defined, 
mesoscale simulations may be performed and structural 
properties calculated for any system of interest. 



B. Mesoscale Simulations of Coarse-Grained 
Polymer Mixtures 

In continuing our multi-scaling approach, mesoscale 
simulations (MS) of various binary polymer mixtures 
with x — were performed using the effective poten- 
tial calculated in the previous section. These simulations 
provide a reference system representative of the case in 
which the blend is well mixed and fluctuations are small. 
Systems investigated were blends of polyethylene (PE), 
polyisobutylene (PIB), and polypropylcncs in their head- 
to- head (hhPP), isotactic (iPP), and syndiotactic (sPP) 
forms. Table [I] shows the simulation parameters used in 
both UA and MS descriptions for each polyolefin blend 
studied. 



TABLE I: Simulation Parameters for Polyolefin Blends (N A = 
N B = 96) 

Blend [A/B] <f) p [sites/ A] R gA [A] R gB [A] 7 



hhPP/sPP 


0.50 


0.0332 


12.18 


13.87 


1.14 


hhPP/PE 


0.50 


0.0332 


12.32 


16.48 


1.34 


PIB/PE 


0.50 


0.0343 


9.76 


16.38 


1.68 


PIB/sPP 


0.50 


0.0343 


9.76 


13.78 


1.41 


iPP/PE 


0.75 


0.0328 


11.33 


16.69 


1.48 


hhPP/PIB 


0.50 


0.0343 


12.41 


9.69 


1.28 



Our MS MD simulations were performed on systems 
of point particles evolving in the microcanonical (N, V, 
E), ensemble. Initially all particles were placed on a cu- 
bic lattice with periodic boundary conditions, where the 
type of particle (A/B) occupying a particular lattice site 
was determined at random. The potential and its cor- 
responding derivative were entered as tabulated inputs 
to the program, and linear interpolation was used to de- 
termine function values not supplied as the algorithm 
proceeds. Each site was given an initial velocity pooled 
from a Mersenne Twister random number generator. |32j 
and the system subsequently was evolved using a veloc- 
ity Verlet integrator. We used reduced units such that all 
the units of length were scaled by R 9 ab (r* = r/R gA B) 
and energies were scaled by k B T. 

Equilibrium was induced by a quenching procedure, 
in which velocities were rescaled at regular intervals to 
maintain the desired temperature. Proper equilibration 
was verified by observing a Maxwell-Boltzmann distri- 
bution of velocities, a steady temperature, a stabilized 
Boltzmann H-theorem function, and a decayed transla- 
tional order parameter. Once velocity rescaling was dis- 
continued during the production stage, we monitored the 
fluctuating steady temperature to ensure the system re- 
mains in equilibrium throughout the duration of the sim- 
ulation. 

During the production stage, trajectories were col- 
lected over a traversal of ~ 8R g . A typical MS MD 
simulation included ~ 4000 particles, evolving for 30,000 
computational steps. Several simulations were run using 
the LONI TeraGrid system [33] to facilitate performing 
numerous simulations at a time. As a benchmark, a typ- 
ical run on a single CPU workstation took ~ 5 h of wall 
clock time; however, since our codes are not yet subject 
to any parallelization process, we stress that this repre- 
sents an underestimate of the potential gain in compu- 
tational time as opposed to running full atomistic or UA 
simulations. 

TABLE II: Mesoscale Simulation (MS-MD) Particle Number 
and Box Dimension Compared to UA Box Dimension. All 
UA simulations are for 1600 chains. 

Blend [A/B] Particles L M s [A] L UA [A] 



hhPP/sPP 


2744 


199.07 


166.61 


hhPP/PE 


5324 


246.21 


166.61 


PIB/PE 


4096 


218.66 


164.91 


PIB/sPP 


5488 


245.19 


164.91 


iPP/PE 


1728 


168.57 


167.27 


hhPP/PIB 


3456 


230.73 


164.91 



Mesoscale simulations provide the center-of-mass to- 
tal correlation functions that describe the polymer mix- 
tures on the large scale and are readily calculated from 
the simulation coordinates. As an exmple we show in 
Figure 1(a) the plot of h AA (k) for a 50:50 mixture of 



6 



hhPP/sPP (\ = 0). Data from mesoscale simulations 
and theoretical predictions are compared against united 
atom simulations for the center-of-mass total correlation 
functions and show an excellent agreement. Analytical 
theory, mesoscale simulations, and united atom simula- 
tions are all consistent in depicting the structure of the 
fluid on the length scale of the polymer radius of gyration 
and larger. Although all of the structural information is 
already contained in the analytic expression, the MS-MD 
simulations are useful as they can provide further infor- 
mation, for example, on the dynamics of the system close 
or far from equilibrium, both of which could be in prin- 
ciple investigated with our coarse-grained systems. 



i 1 i 1 r 
(b) 




0.1 0.2 0.3 0.4 



k(A ) 



0.1 0.15 0.2 0.25 

k (A 1 ) 



FIG. 1: (a) Plot of h% A (k) for hhPP/sPP obtained from 
mesoscale simulation (open red circles). Comparison with 
UA data (filled circles) and theoretical predictions (solid line) 
shows quantitative agreement, (b) Plot of h™™(k) calculated 
using the inverse mapping procedure, Equation 11 (open red 
circles) compared to data from the full UA MD simulation 
(solid circle). 



However, while h cc (k) describes the center of mass cor- 
relations between particles, when we are concerned with 
the structure of the liquid on the local scale, it is nec- 
essary to look at the monomer total correlation func- 
tions. Although the large scale information is completely 
provided by the mesoscale information, the information 
on the local-scale structure is averaged out during the 
coarse-graining. To take advantage of the gain in com- 
putational time given by mesoscale simulation, it is con- 
venient to combine the large scale information from the 
latter with a short united atom simulation that describes 
the local structure. The two are merged together through 
a multiscale modeling procedure as presented in the next 
section. The advantage in performing this two-steps pro- 
cedure is that the simulation with atomistic resolution 
can now be limited to short simulation runs and to a 
small ensemble of molecules, as the needed large-scale, 
long-time information comes from the fast mesoscale sim- 
ulations. 



C. Recovering the Monomer Level Description: 
Implementation of the Multiscale Modeling 
Procedure for Polymer Blends 

In our multiscale modeling we combine the mesoscale 
simulation, which correctly captures the large scale, 
global structure, with UA simulations, needed to deter- 
mine the local properties of the polymer. By this pro- 
cedure we aim at bridging the different length scales 
characteristic of the mixture and obtaining a complete 
description of the polymer blend structure at a defined 
temperature and composition. 

The key to this procedure lies in extracting the 
monomer-monomer correlation function, h mm (k) from 
the coarse-grained correlation function, h cc (k), by inverse 
mapping, using the reciprocal relation: 



'<q m (fcvgr(*o 



(ii) 



Here, h cc {k) is obtained from the mesoscale simulation, 
while the intramolecular form factors, w(fc), entering 
Equation II are calculated directly from short UA MD 
simulations. Since UA-MD simulations are needed only 
to obtain the intramolecular distributions, these simu- 
lations need only be for system box sizes on the order 
of the radius of gyration, which is a significant advan- 
tage over obtaining the global structure from UA-MD 
alone which requires much larger box sizes to accurately 
calculate the intermolecular distributions. The calcu- 
lated h mm {k) is valid only for small k values, where the 
coarse-grained description applies, and begins to diverge 
as u cm (k) approaches zero at large k. This is shown in 
Figure 1(b), where we compare directly h mm (k) calcu- 
lated using Equation [IT] with UA-MD data of the same 
system. 

To include local scale information, the total correla- 
tion function h mm (k) at small k, obtained following the 
described procedure, is combined with h mm (k) from UA 
simulation for large values of k. The key point in the 
procedure is to define a cut-off length, k cu t oc 1/ r cut , at 
which the two descriptions are to be combined. 

Because the local scale information is averaged out in 
the coarse-grained description, it is important to ensure 
that the whole intramolecular description is accounted 
for by the UA simulation. To estimate the length scale 
at which local intramolecular contributions become neg- 
ligible, we evaluate the fraction of intra to total site/site 
contacts. For a given component of type a in the mix- 
ture, this function is defined as 



fsaa (0 



N saa (r) 
N tota i(r) 



(12) 



where N saa (r) is the number of a type intramolecular 
contact sites 

N saa (r) = 4np f \r') 2 0J™™{r') dr' , (13) 
Jo 
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while the total number of site/site contacts is given by 



III. RESULTS 



Ntotaiir) = -npr 3 



(14) 



In our procedure, the integral in Equation [13] is computed 
numerically using UA data for w™™(r). 

Selecting the correct cut-off distance is an important 
step of the procedure. Choosing a larger cutoff radius 
r cu ti smaller k cutl includes more information from UA 
simulations, increasing the precision of the calculation, 
but it increases the computational time, partially de- 
feating the purpose of adopting a multiscale procedure. 
In the case of polymer melts we observed that a value 
°f fs( r ) = 0.025 gives good precision in the calculated 
h mm (r) when compared with full united atom simula- 
tion data, while retaining a reasonable speeding up of the 
calculations. [18] For polymer blends we adopt the same 
criterium and we evaluate a posteriori the agreement ob- 
tained for the blend samples analyzed in this paper. In 



Table III we report, for both AA and BB type inter- 
actions, the effective radius, r cut , for combining simula- 
tions, which lies in the intermediate length scale on the 
order of a few R g units and is unique for each blend. Also 
reported is the total number of site/site contacts within 
the radius determined by the value f s (r) = 0.025, evalu- 
ated using Equation [14] at r = r cut . Dividing N tota i by 
the number of united atoms per chain (96 in our case) 
gives the number of chains, n, needed in a UA simulation 
to produce the statistical information necessary to obtain 
the local structure. This is important since it determines 
how "short" the UA-MD simulations must be without 
losing pertinent information about local correlations. 

Large-scale and local-scale information were combined 
at the chosen k cut = 27r/r cut , and the procedure was re- 
peated for each of the mixtures shown in Table [I] The 
correlation coefficient between h mm (r) determined from 
our multiscaling procedure and from full UA simulations 
is calculated for the case when the multiscale simulations 
are combined at f s = 0.025. These values are also pre- 
sented in Table [Till showing that once r cut is defined, the 
multiscale procedure provides an accurate way of obtain- 
ing the total correlation function for binary blends while 
circumventing the need for prohibitively large atomistic 
simulations. For the case of AB type interactions, an av- 
erage radius between AA and BB data is used to deter- 
mine where the simulations should be combined. While 
this approach allows one to explore a large range in the 
degree of polymerization, the method becomes impracti- 
cal for liquids of long chains, for which intramolecular ef- 
fects will remain long-ranged, since the cutoff length scale 
occurs on the order of R g . Moreover, for long chains, en- 
tanglement effects enter the dynamics and have to be ac- 
counted for. Therefore, for large N systems it is advisable 
to include a third, intermediate, level of coarse-graining, 
with a cut-off of the order of (and larger than) the per- 
sistence length and the entanglement length scale. 



A. Total Pair Correlation Functions of Athermal 
Reference Blends 

The method just discussed gives a systematic way of 
merging simulations to optimize the tradeoff between the 
gain in accuracy due to inclusion of UA simulation data 
and the gain in efficiency due to the coarse grained meso- 
scopic picture. This procedure works well as it yields to- 
tal correlation functions in excellent agreement with UA 
simulations at a reduced computational cost. 



TABLE III: Effective Radius for Combining Simulations and 
the Number of Total Site/Site Contacts Evaluated at f s (r) = 
0.025 



System Type r c „ t L4]/ 8 



N, 



tot a I 



n Corr. Coeff 



hhPP/sPP 
hhPP/sPP 
hhPP/PE96 
hhPP/PE96 
PIB/PE96 
PIB/PE96 
PIB/sPP 
PIB/sPP 
iPP/PE96 
iPP/PE96 



AA 
BB 
AA 
BB 
AA 
BB 
AA 
BB 
AA 
BB 



29.1 
28.7 
28.7 
27.6 
28.8 
27.2 
29.1 
29.6 
29.3 
27.6 



3427 
3288 
3288 
2924 
3432 
2891 
3536 
3721 
3458 
2891 



36 
34 
34 
31 
36 
30 
37 
39 
36 
30 



0.9999 
0.9999 
0.9998 
0.9999 
0.9999 
0.9999 
0.9995 
0.9999 
0.9999 
0.9998 



Results from the multiscale procedure are reported in 
Figures [2] - [6] In Figure [2j the left panels show the to- 
tal correlation function, h mm (k), obtained by combining 
UA simulations for the short length scales (large k) with 
mesoscale simulations for the long length scales (small k) , 
for the first mixture in Table [TJ hhPP/sPP. The vertical 
dashed line in the left panels represents the cut-off length 
scale at which the two simulations are combined, which 
corresponds to f s = 0.025. The insert in the left panels 
depicts the peak representing the local structure, which 
is accurately determined by UA simulation. This peak 
depends on the geometry of the monomeric structure and 
changes as a function of the type of polymers involved in 
the mixture. 

The right panels in Figure [2] show h mm (r), i.e. the 
Fourier transform of the corresponding total correlation 
function, h mm (k). As in the case for melts, we adopted a 
sampling step of Afc = 0.01 in reciprocal space, as this is 
of the same order as the discontinuity in h mm (k), which 
results from joining the two simulations. Because of the 
amplitude of the chosen sampling step, there is no ef- 
fect on the Fourier transformed function h mm (r). The 
pair distribution function calculated from the multiscale 
modeling procedure is identical to the coarse-grained an- 
alytical expression in the large-scale regime, but also in- 
cludes the local scale solvation shells, which come from 
UA-MD and are not captured by the PRISM thread ex- 
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pressions. Both pieces of information are needed to pro- 
vide the complete structural and thermodynamic picture 
of the system. 

Figure [2] illustrates the spirit of our multi-scale model- 
ing approach in which independent simulations represent- 
ing the same system at different levels of coarse-graining 
can be combined to provide a complete description of 
the polymer. Analogous plots for the hhPP/PE and the 
PIB/PE mixture are reported in Figure [3] and Figure 
|4j respectively. Finally, results for the two mixtures of 
PIB/sPP and iPP/PE are presented in Figures [5] and [6] 
All systems in real space show quantitative agreement 
with UA-MD data but are obtained at a much more 
efficient computational time than running the full UA 
MD simulation. Furthermore, the procedure to obtain 
the pair correlation functions is entirely analytical, hence 
we do not utilize any optimization procedure or numeri- 
cal fitting scheme to obtain consistency between the two 
descriptions, and thus the method is fully transferable. 
Figures [2][6] demonstrate the versatility of the approach 
for studying mixtures of polymers with subtly different 
chemical architecture, since the same multiscale proce- 
dure is applied in all cases. 



B. Applications to Thermal Blends with Realistic 
X Parameters 

In the multiscale modeling procedure presented above, 
we set x = in ° ur soft-colloidal representation, and 
the UA simulations to which we combine our mesoscale 
simulations are assumed to be well-mixed. This provides 
an efficient means of obtaining pair correlation functions 
for polymer mixtures close to the athermal limit. In this 
section we examine the assumption that the blends are 
well-mixed and discuss the application of our multiscale 
modeling procedure to modeling thermal mixtures where 
subtle local features such as monomer shape, branching, 
and site energetic asymmetries gives rise to non-trivial 
divergent fluctuations. 

Blends of polypropylenes of different tacticities have 
been studied by Woo et al., and the x parameter of blends 
of iPP/sPP and aPP/sPP was found to be nearly zero 
over the temperature range T=423-453 K.[33] Thus, for 
the hhPP/sPP blend of Figure [2] it is reasonable to as- 
sume the effective x parameter will be small, which is 
supported by the overall quantitative agreement in Fig- 
ure [T] and [3j For the other systems investigated in this 
publication, many of the pertinent x parameters may 
be determined from the available literature on polyolefin 
blends which is presented in Table |IV| 

To demonstrate the extension of our approach to ther- 
mal mixtures, we use the literature values for the x P a ~ 
rameter evaluated at the simulation temperature, T = 
453K and perform MS-MD simulations for the mixtures: 
hhPP/PE, PIB/PE, and hhPP/PIB. The hhPP/PIB 
blend is particularly interesting because it exhibits a 
lower critical solution temperature (LCST) phase dia- 



k(A') 







k(A') 



r(A) 



FIG. 2: (Left) Multiscale modeling: The left panels show the 
total correlation function, h mm (k), for A A (top), BB (mid- 
dle), and AB (bottom) interactions, for a mixture of 50:50 
hhPP/sPP. The data over the range of small k values was ob- 
tained by mesoscale simulation, whereas over the large k range 
it was obtained by UA MD simulation. The inset depicts the 
local structure. The dashed line indicates the value at which 
the two simulations were combined. (Right) The correlation 
function, h mm (r), after Fourier transform (solid red line) is 
compared with results from the full UA MD simulation (open 
symbols) . 



gram, which can in principle be modeled using our mul- 
tiscale approach as well. [IT] Since the parameters used in 
our mutiscale modeling procedure are determined using 
a UA description where each CH^ group is represented 
as a effective site, care must be taken when calculating 



the effective parameter from Table IV where x is de- 
fined on a monomer basis, for which several CH^ groups 
are grouped as one monomer unit. Since x is propor- 
tional to the free energy of mixing per site, the x P a ~ 



TABLE IV: Temperature Dependence of Polyolefin Blends 



Blend [A/B] 


X(T) 


X/Xs (453^) 


hhPP/PE 


-0.0294+ 17.58/T" 


0.077 


PIB/PE 


0.00257 + 4.99/T 6 


0.163 


iPP/PE 


0.01 c 


0.360 


hhPP/PIB 


0.027- 11. 4/T d 


0.018 



Reference 26 b Reference 27 c Reference 25 d Reference 23 
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FIG. 5: Same as in Figure[2]for a blend of 50:50 PIB/sPP. 
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FIG. 4: Same as in Figure[2]for a blend of 50:50 PIB/PE. FIG. 6: Same as in Figure[2]for a blend of 75:25 iPP/PE. 
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rameter must be normalized by the average number of 
united atom sites per monomer [4] as defined in the lit- 
erature (6 for hhPP/PE 35J, 4 for PIB/PE[3B], and 4.8 
for hhPP/PIB0]). For hhPP/PE this corresponds to a 
value of x = 0\0016; for PIB/PE a value of x = 0.0034; 
and hhPP/PIB a value of 3.8 x 10~ 4 . For the temper- 
ature range at which the united atom simulations were 
performed, the magnitude of x m an cases is small, sup- 
porting our initial modeling of \ — in the previous 
section. In Table |IV| we present the thermodynamically 
relevant parameter, x/Xs at 453-ftT, which provides an 
indication of how far the mixture is from the spinodal. 

We calculated the total correlation function for 
the thermal mixtures of hhPP/PE, PIB/PE, and 
hhPP/PIB, using the multiscale modeling approach de- 
scribed above. These results are presented in the sup- 
plemental material|37] and look nearly identical to the 
athermal results shown in Figures [| and g] Although 
the correlation functions are nearly identical with the 
athermal limit because the blends are well mixed, sub- 
tle differences in the structure as a result of increased 
concentration fluctuations should be manifest in the dis- 
tributions. To assess the changes in structure that re- 
sult from increased fluctuation, the monomer level par- 
tial static structure factors can be calculated from the 
total correlation function, h mm (k), 

ST A \k) = ^Z^ + ptfhJZik) , 

S£b (*0 = (1 " <t>)"BB (*0 + p(i - <l>) 2 h<EZ(k) ,(15) 

s%E l (k) = P <t>{i-cj>)hj%{k) , 

where the monomer form factors were determined from 
UA simulations as in Equation [IT] above. The structure 
factor measuring correlations in the relative concentra- 
tion, S^{k), which diverges as the mixture approaches 
the spinodal, is expressed as a linear combination of these 
partial structure factors, 



30 - 



fS^{k)+<i> 2 S^{k)- 



■24>{l-<t>)S^{k). 



S**(k) = (1- 

In small angle neutron scattering (SANS) experiments, 
the x parameter is determined from fitting the par- 
tial structure factor to the random phase approximation 
(RPA) equation of de Gennesf3"5] 



1 



+ 



s(k) <i>uyz{k) (i-0) w ™™(fc) 



2X, (17) 



where for convenience the monomer site volumes were set 
equal to one. 

Figure [7] presents S^(k) for the three different ther- 
mal blends obtained from our multiscale modeling pro- 
cedure. The static structure factor calculated in this 
manner exhibits good agreement when compared to the 
RPA equation, which was evaluated at x — 0.0016 for 
hhPP/PE, x = 0.0034 for PIB/PE, and x = 3.8 x 10~ 4 
for hhPP/PIB, using the intramolecular form factors 
from UA MD simulations. The relevance of this com- 
parison to the RPA formula, which is known to fit low 



A) 



X = 0.0016 

x = 0.00 




tt S88 8 



k(A ) 



30 


- B) 


1 1 1 1 1 
\ — x = 0.0034 


25 




*\ — x = 0.00 


20 


~*\ \ 

" \ D 
\ \ 

- M 

\ V) 






\ \ 

\\ 


\ - 


t 15 


\\ 

- \\ 
\\ 

\\o 

\\ 
\\ 




10 

5 


- <P 
\ 


V 







1 , 1 1 1 





0.2 0.4 0.6 






k (A"') 


30 




1 1 1 1 

— X = 0.00038 - 


25 




— x = o.oo 




- *k 




20 






2~ 






t 15 




•\ 
•\ 


10 






5 











11 1***8888 

i , 1 , 1 



0.6 



k(A ) 



FIG. 7: (A) The concentration fluctuation structure factor 
obtained from the multiscale modeling procedure (open cir- 
cles) for hhPP/PE (<f> = 0.5) at T = 453K. For compari- 
son the RPA equation evaluated at \ = 0.0016 (solid line) 
and x — 0-00 (dashed line) is shown. Filled circles represent 
the structure factor from the full UA simulation. (B) The 
same as part (A), except for PIB/PE (<j> — 0.5), for which 
the RPA equation, evaluated at \ — 0.0034 (solid line) and 
X = 0.00 (dashed line), is shown. (C) the same except for 
the hhPP/PIB blend, and the RPA equation is evaluated at 
X = 0.00038 (solid line) and x = 0.00 (dashed line). 



wavevector scattering curves well, is that there is clearly 
a better fit using the experimental x parameter than with 
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the x = case, which demonstrates the consistency of the 
description and that our approach is able to capture the 
fluctuations in concentration that arise in thermal poly- 
mer mixtures even at the relatively high temperatures of 
these simulations. Also, the advantage of a multiscale ap- 
proach is exhibited by Figure [7] since the low wave vector 
regime is determined by mesoscale simulations so that 
UA simulations only need to be performed on systems at 
length scales up to the cut-off length (dashed line) . This 
is important since only the initial stages of the divergent 
behavior need to be captured by united atom simulations, 
thus the need for prohibitively large simulation boxes is 
circumvented. For comparison, S(k), calculated for the 
full UA MD simulation is also shown (offset for clarity) 
and agrees with our multiscale results, demonstrating the 
consistency between the two approaches. 

When experimental values of the x parameter are un- 
available, a commonly used alternative is to estimate the 
interaction parameter from the solubility parameters cal- 
culated from a pure melt. These may be obtained from 
the cohesive energy determined from MD simulations of 
the individual melts. For the systems investigated here, 
we used united atom simulations of pure melts provided 
by Jaramillo et al.0] The cohesive energy density can be 
calculated from the radial distribution functions for each 
UA site type, g^, where i, j designate a particular C, CH, 
CH 2 , or CH 3 group, according to 



parameter may be estimated according to 

(5 A - 5 B ? 



XH 



(19) 



where the densities, p° , are the pure melt densities of 
species A or B. Table [Vl| shows a comparison between the 
Flory-Huggins x parameter determined from this conven- 
tional solubility parameters approach and the observed 
experimental values. A major limitation of this solubil- 
ity approach is that the value of x is always positive 
and UCST phase behavior is always predicted. Further- 
more, because non-combinatorial entropy is ignored, the 
concentration dependence of x cannot be predicted. De- 
spite these assumptions, the solubility approach has been 
demonstrated to hold up reasonably well for a large vari- 
ety of polyolefin blends |H]. For comparison, in Table VI 
we also present the x parameter computed directly from 
the full UA simulation of a polymer blend, by an extrap- 
olation of the three lowest wave vector points to k = on 
an Ornstein-Zernike plot (1/S(k) vs. fc 2 ), analogous to 
that obtained in typical SANS experiments. The lowest k 
vector is determined by the size of the simulation box and 
in this case is k — 0.4 A. Due to this limitation, the effec- 
tive x parameter can be shifted by about ±0.002 while 
still maintaining a reasonable fit to the blend structure 
factor. [5] 



U, 



coh 



2tt 



J2 prfj 

ij 



u ij( r )9ij( r ) r2 dr, (18) 



where Vij is the potential between nonbonded sites, for 
which we employ a Lennard- Jones potential with the 
corresponding TraPPE-UA parameters of Martin and 
Sicpmann|39j. In Equation 18 the integration is carried 



out over the attractive branch of the potential. The sol- 
ubility parameter (5 = \J — C/ co h), calculated for each of 
the polyolefins in this study is presented in Table |V| For 
comparison we also present results obtained by Piitz, et 
al.[40j from UA simulations of chains with 24 backbone 
carbons. 



TABLE V: Solubility Parameter for UA simulations of poly- 
olefin melts (N=96) at T — 453K. Values for Smd-2 are 
computed from Ref. [40, for chains with 24 backbone car- 
bons. 

Polymer 8 (MPa)2 5 MD - 2 



PE 


14.4 


13.7 (T = 


hhPP 


13.4 


13.3 


iPP 


12.5 


12.4 


PIB 


13.9 


14.1 


sPP 


13.0 


12.2 



TABLE VI: Comparison between experiment, solubility pa- 
rameter, and RPA approaches to determine the \ parameter 
(T = 453K) for the polyolefin blends in Figure [FJ 



Blend [A/B] 




XH 


Xoz plot 


hhPP/PE 


0.0016 


0.0048 


0.0022 


PIB/PE 


0.0034 


0.0012 


0.0027 


hhPP/PIB 


0.00038 


0.0012 


0.0017 



Finally we note that much work has been done using 
the PRISM theory of Schweizer and Curro to determine 
both the qualitative and quantitative miscibility behav- 
ior of polymer blends. [42 44J For polyolefin blends of the 
type investigated here, realistic models utilizing PRISM 
expressions for the pair distribution function along with 
UA potential parameters have been shown to give rea- 
sonable estimates of x f° r various polypropylene and 
polyethylene blends. |5} |45[ [46] ; however this approach is 
not pursued further here since the purpose of the current 
publication is to demonstrate the implementation of the 
multiscale modeling procedure once the x parameter is 
obtained. 



IV. CONCLUSION AND OUTLOOK 



From these solubility parameters, the Flory-Huggins 



In this paper we present a multiscale modeling pro- 
cedure to simulate mixtures of polymer chains at large 
scales. The method is completely general and applies to 
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mixtures of polymers with different molecular structures, 
at different thermodynamic conditions of temperature, 
mixture composition, and chain length. The procedure 
combines large-scale information from a mesoscale sim- 
ulation with a more detailed model which captures the 
local structure of the blend. For the detailed model, we 
use a united atom representation which captures the lo- 
cal polymer structure. Our mesoscale model represents 
the blend as a mixture of soft colloidal particles and 
accurately describes the large-scale structure as exhib- 
ited in the center of mass radial distribution function. 
Once the mesoscale simulations are performed, the rele- 
vant chemical details are reinserted by implementation of 
the Ornstein-Zcrnikc formalism. We then combine this 
information with short united atom simulations to ob- 
tain a complete description over all length scales of in- 
terest. In this article we test the approach by applying 
our method to several different polyolefin mixtures. We 
check the validity of our procedure by direct compari- 
son of the pair correlation function against full united 
atom molecular dynamic simulations. Because the UA- 
MD simulations available for comparison are limited to 
high temperatures where the samples can be more eas- 
ily equilibrated, the test of our procedure includes only 
high temperature samples. However our mesoscale simu- 
lations are easily performed at temperatures approaching 
the demixing temperature. Furthermore, the multiscale 
modeling procedure provides a straightforward method 
to circumvent large simulations when modeling thermal 
mixtures, which has the potential to extend the range 



currently available to computer simulations of polymers. 
At low temperatures, where equilibration of even the lo- 
cal UA-MD simulation may be difficult, our two step pro- 
cedure could be implemented "on the fly" where informa- 
tion is transferred from the local UA-MD (which provides 
the packing information) to the global scale MS-MD and 
vice versa until consistency between the two is reached. 

One advantage of the proposed method is that no op- 
timized parameterization scheme is required to maintain 
consistency between the different levels of description. 
The approach is thus transferable to different systems 
and different conditions, while remaining quantitative 
enough to distinguish even subtle differences in chemical 
structure. Considering the advantages of the proposed 
multiscale approach, it should provide an important tool 
for investigating the dynamics of large-scale phenom- 
ena via MD simulation. Future work should focus on 
performing simulations of large systems approaching the 
spinodal where full UA MD simulations are prohibitively 
difficult. 
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